Formulations of the 3+1 evolution equations in curvilinear coordinates 
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Following Brown in this paper we give an overview of how to modify standard hyperbolic 
formulations of the 3+1 evolution equations of General Relativity in such a way that all auxiliary 
quantities are true tensors, thus allowing for these formulations to be used with curvilinear sets of 
coordinates such as spherical or cylindrical coordinates. After considering the general case for both 
the Nagy-Ortiz-Reula (NOR) and the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulations, 
we specialize to the case of spherical symmetry and also discuss the issue of regularity at the origin. 
Finally, we show some numerical examples of the modified BSSN formulation at work in spherical 
symmetry. 
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I. INTRODUCTION 

In 3+1 formalism of General Relativity one splits 
spacetime into a foliation of 3-dimensional (3D) space- 
like hypersurfaces (assuming that the spacetime is glob- 
ally hyperbolic) , and projects the Einstein field equations 
in the normal and tangential direction to those. In this 
way, the 10 independent field equations are naturally sep- 
arated into 4 constraint equations and 6 evolution equa- 
tions for the geometric degrees of freedom. The evolution 
equations that are obtained directly from this projection 
are known as the Arnowitt-Deser-Misner (ADM) equa- 
tions dfl. As was already realized in the late 80's and 
early 90's, these ADM evolution equations, though phys- 
ically correct, have nevertheless one serious drawback: 
they turn out to be only weakly hyperbolic and as such 
are not mathematically well-posed (see e.g. 0). By this 
one means that the solutions do not depend continuously 
on the initial data and can be unstable in the presence 
of constraint violations, which in practice implies that 
one will encounter serious stability problems in numeri- 
cal evolutions based on these equations. 

It turns out, however, that one can construct alterna- 
tive formulations of the evolution equations by adding to 
them multiples of the constraints in a variety of differ- 
ent ways. These new systems of evolution equations will 
have the same physical (constraint satisfying) solutions, 
but will typically differ significantly in their mathemat- 
ical structure. Over the last two decades, a number of 
different well-posed strongly hyperbolic formulations of 
the 3+1 evolution equations have been proposed, and 
several of them have been tested in numerical evolution 
codes [2i| . In particular, the formulation proposed by 
Shibata and Nakamura, and Baumgarte and Shapiro, 
known as the BSSN formulation [f| |(|, has turned out 
to be very stable and robust in practice, and has become 
the standard formulation used by most 3+1 evolution 
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codes today. This formulation has finally allowed the 
accurate simulation of binary black-hole systems with 
different masses and spins, starting from wide separa- 
tions through the merger and ring-down of the final black 
hole 0-0 (one should mention, however, that the first 
successful simulation of multiple orbits of binary black 
holes was in fact carried out by F. Pretorius using a very 
different approach based on the so-called generalized har- 
monic formulation, an approach that is still being used 
today by a number of different groups [irj||). 

The BSSN formulation, though very successful in prac- 
tice, has the drawback of involving dynamical quanti- 
ties that are not true tensors, such as tensor densities 
and contracted Christoffel symbols. This represents no 
problem in most 3D simulations where one typically uses 
Cartesian coordinates, but becomes an important issue 
when one considers curvilinear coordinate systems, such 
as spherical or cylindrical coordinates. 

Recently, Brown introduced a more general version of 
the BSSN system where all dynamical quantities are true 
tensors [l|. This "generalized BSSN" formulation is thus 
ideally suited for the use of curvilinear systems of co- 
ordinates, which in particular allows one to construct a 
BSSN version of the evolution equations for the case of 
spherical or cylindrical symmetry. 

In this paper we give an overview of the main ideas 
behind Brown's approach, and apply them to both the 
Nagy-Ortiz-Reula (NOR) (nj and BSSN formulations. 
The paper is organized as follows. In Section [TT| we give 
a brief review of the 3+1 formalism. Later, in Section Hill 
we discuss some important results related to the fully co- 
variant expressions of the Riemann and Ricci curvature 
tensors in terms of a background metric. Section ITVl then 
considers the case of the NOR formulation and its gen- 
eralization to curvilinear coordinates. In Section [V] we 
repeat the same analysis for the BSSN formulation, and 
also include a brief discussion of the Gamma driver shift 
condition. In Section IVIl we consider the particular case 
of BSSN in spherical symmetry, and discuss the basic 
equations and the important issue of the regularization 
at the origin. Finally, in Section IVIII we present some 
numerical examples. We conclude in Section IVlIII 
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Throughout the paper we will use geometric units such 
that G = c = 1. Also, Greek indices will represent all 
spacetime dimensions and will run from to 3, while 
Latin indices will represent only spatial dimensions and 
will run from 1 to 3. 



II. BASIC 3+1 EQUATIONS 

Before considering the NOR and BSSN formulations 
it is convenient to first review the basic concepts and 
equations of the 3+1 formalism of general relativity (for 
a more detailed introduction see e.g. [2]). 

In the 3+1 formulation spacetime is foliated into spa- 
tial hypersurfaces parametrized by a time function t. 
The basic dynamical quantities of the 3+1 formulation 
are then the metric of the spatial hypersurfaces 7y and 
the extrinsic curvature tensor of those hypersurfaces 
which is defined as 

:= -i£ V a n y (2.1) 

where is the time-like normal vector to the spatial 
hypersurfaces, and :— Sp + n a ng the projection op- 
erator onto the hypersurfaces. 

Furthermore, one also introduces the lapse function a 
that measures the proper time elapsed between adjacent 
hypersurfaces along the normal direction, and the shift 
vector ft 1 that controls how the spatial coordinates prop- 
agate from one hypersurface to the next. In more detail, 
an observer moving along the normal direction to the 
hypersurfaces (also known as an Eulerian observer) will 
have a coordinate speed given by — /3 l , and will measure 
a proper time dr = adt. In terms of these coordinates, 
the unit normal vector becomes = (1/a, — (3 l /a), and 
the extrinsic curvature tensor takes the form 

Kij = (drfij - £p7ij) , (2.2) 

with £p the Lie derivative with respect to the shift vec- 
tor, and where we have only considered spatial compo- 
nents using the fact that the extrinsic curvature is by 
definition normal to the hypersurfaces. 

Given the spacetime foliation just described, the Ein- 
stein field equations separate naturally into two distinct 
groups. The first group corresponds to those equations 
that have no time derivatives and results in the so-called 
Hamiltonian and momentum constraints 

H := - (R + K 2 - l\,,l\' J ) - 8tt P = , (2.3) 
M i := Vj (K ij - 7 lJ K) - 8nf = . (2.4) 

In the above equations R := 7 iJ R^ is the trace of the 
spatial Ricci tensor R^, K :— j 1 ^ Kij is the trace of the 
extrinsic curvature, and V, is the covariant derivative 
associated with the spatial metric 7y, while p and j l 
are the energy and momentum densities measured by the 



Eulerian observers and are given by 

p := nWT^ (2.5) 
f := -P^n v T^ , (2.6) 

where T^ v is the stress-energy tensor of the matter. 

The second group of field equations corresponds to 
the true evolution equations of the system. In terms of 
the quantities introduced above these evolution equations 
take the form 

dtlij - £$1ij = -2aKij , (2.7) 
d t Kij - fpKij = -VjVj-a 

+ a ^ + KKij - 2K ik K)\ 

+ 47ra [jij (S-p)- 2Sij] , (2.8) 

where Rij is the 3-dimensional Ricci tensor associated 
with the spatial metric 7^ : 

Rij = — ^ l^dmdnlij + lfm(idj)T m + r"T(y) m 

I 9pnm p _i_ p pmn /n q\ 

"t" zl (2 1 j)mn 1 1 mni L j ) ) 

with r* := 7 mn r l m „, and where Sij is the stress tensor 
measured by the Eulerian observers defined as 

Sij := P?P?T a p , (2.10) 

with 5* := j^Sij. The evolution equations above are 
known in the numerical relativity community as the 
Arnowitt-Deser-Misner (ADM) equations [1, |j|- 



III. CURVATURE TENSOR IN TERMS OF A 
BACKGROUND METRIC 

It is convenient at this point to review some well-known 
fully covariant expressions for the Riemann and Ricci cur- 
vature tensors in terms of a background metric. 

Let us assume that we have a manifold with some co- 
ordinate system and two different metric tensors defined 
on it: the "physical" metric 7^, and some "background" 
metric %j that is not necessarily flat (though in the fol- 
lowing sections we will assume that the background met- 
ric is indeed flat). We now want to express the curvature 
tensor associated with the physical metric 7^ in terms of 
the curvature associated to the background metric 7^, 
together with covariant derivatives of 7y with respect to 
this background. In order to do this, we start by defining 
the quantity: 

A a t c ■— r a fc c — T a bc , (3.1) 

with r a b c and f a b c the Christoffel symbols associated 
with jij and 7^ respectively. Notice that even though 
neither T a i, c nor T a b c are components of tensors, their 
difference A a h c is in fact a proper tensor. 
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Having denned A% c let us now calculate the covariant 
derivative of the physical metric 7^ in the background 
geometry, that is V a 7b c . Notice first that, in general 



V a 7b c = V Q 7 bc = , 

y a 7bc 0. 



(3.2) 
(3.3) 



If we now take the convention that indices of A a (, c are 
raised and lowered with the physical metric 7^ , then we 
can use (13.21) to show that: 



Va7bc = 2A( bc 



and equivalently 



V a7 



be 



-2A {bc K 



(3.4) 



(3.5) 



One can now solve for A a oc from the above expressions 
to find 



A a 6c = - 7 am (Wbjcm + V c 76m - V m7&t 



(3.6) 



Notice that this expression for A a b c is in fact identical to 
that for the Christoffel symbols T a b c , but with the par- 
tial derivatives replaced with covariant derivatives on the 
background. In particular, if the background is flat and 
we use Cartesian coordinates we will have A a f, c = T a bc 
and V a = d a , so that the last expression reduces to the 
standard definition of the Christoffel symbols. 

We can now use (|3.6[) to show that the physical Rie- 
mann curvature tensor can be written in terms of the 
A a bc as: 



R a bcd = R a bcd + 2V[ c A a d ]f, + 2A a TO [ c A r 



(3.7) 



where R a bc d is the curvature tensor of the background. 
Again, the second term has the same structure as the 
standard expression for the Riemann tensor, but with 
the r a h c replaced with A a h c , and the partial derivatives 
replaced with covariant derivatives on the background. 
The expression again reduces to the usual one for a flat 
background in Cartesian coordinates. 

Next, let us lower the first index in the Riemann ten- 
sor. This is not as trivial as it might seem since now 
7y can not be brought inside the operator V a directly. 
After a somewhat lengthy algebra, where one needs to 
use the expression for the commutator of the covariant 
derivatives of a rank 2 tensor in terms of the Riemann, 
one finally finds that: 



bed 



lam 1 *' bed ibrn 1 *' acd 
ac 

V (i V a 7fc c - V c V a 7 6d 

A Am _ a Am 
mad<-^ be ^raac^ bd 



(3.8) 



Notice that in the last expression we do not lower the first 
index of R a bcd with 7^, since by convention it should be 
lowered with 7^- . 



Finally, let us find the expression for the Ricci tensor 
Rab '■= l cd Racbd- Using (I3.8p we find, after some algebra: 

Rab = -^7 m ™ V OT V n 7„ 6 + V„V 6 7 m „ 
- V a V m 7b n - VhV r „7 a 

_i_ A \mn A Am 

' <-±rnna L -* b L -*mab L -* 

-i mn R C mn(a lb)e , (3.9) 

where we have defined 

Am . „,fli Am T^m abf^m /o in\ 

■ = 7 A ab = 1 -7 1 ab ■ (3.1U) 

We can in fact rewrite the Ricci tensor in terms of deriva- 
tives of the quantity A m just defined. Using (I3.6[) one 
finds that (|3.Q[) is entirely equivalent to 



R a 



, = -i7" m V m V„7a fc + 7m (a V fc) A m + A m A (afc)r 



A mn a 1 Amn a 

LL ± {a^b)mn + ^* a^-mnb 

mnnc 
I - tt mn(a lb)c ; 



(3.11) 



For a flat background in Cartesian coordinates, this last 
expression clearly reduces to the standard expression 
given in (|2.9p . 



IV. THE NOR FORMULATION 

A. Standard formulation 

The Nagy-Ortiz-Reula (NOR) formulation [ll[ is in 
essence a generalization of the Bona-Masso formulation 
(BM) of the early 1990s [TJGJ]. This formulation is 
based on first writing the three-dimensional Ricci tensor 
that appears in the ADM evolution equations as: 



R lJ = -\i mn d m d nll3 +7fc( i 7 - ) r fc + r fe r (lj)i 

1 (i*- j)mn 1 ^ jl rnnj i 



(4.1) 



where T l := 7 m ™r i 1 „. 

The crucial difference with the ADM formulation is the 
fact that the quantities T l that appear in the Ricci tensor 
above are now promoted to independent quantities and 
evolved separately. To find the evolution equations for 
the r ! we first note that from their definition we have 



(4.2) 



with 7 the determinant of 7^ . From this one can easily 
show, after some algebra, that: 



+ 1 m d m {aK) 



(4.3) 
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where the Lie derivative of T l that appears in the last 
expression should be understood as that of a vector: 

£pP := P m d m T i - r m d m /3* . (4.4) 

In fact, one can now add a multiple of the momentum 
constraints (|2.4[) to this equation to obtain a final evolu- 
tion equation of the form: 

d t T l - £pT l = 7 " m a m 9„/3 l - [2K im - Y m K] V m a 

+ 2 aK mn T\ nn - 87ra£f , (4.5) 

with £ an arbitrary parameter. The importance of adding 
a multiple of the momentum constraints to the evolution 
equation for T* comes from the fact that, if one chooses 
a slicing condition of the Bona-Masso family 

d t a — £pa = —a 2 f(a)K , (4.6) 

with /(a) an arbitrary function of a, then the NOR for- 
mulation can be shown to be strongly hyperbolic (and 
thus well-posed) if one takes £ = 2 and / > 0, or more 
generally if one takes £ > 0, / > and / ^ 1 (see e.g. 
reference Q). 

Instead of using the Bona-Masso slicing condition, 
one can assume that the densitized lapse defined as 
q := a~/~^ 2 , with / a constant, is an a priori known 
function of spacetime, a = F(t,x l ). The same results 
about hyperbolicity then follow. 

The standard NOR formulation in fact also adds an 
arbitrary multiple of the Hamiltonian constraint of the 
form arj^/ijH to the evolution equation of Kij, with r\ 
another free parameter (for rj ^ one finds a new region 
of parameter space where the system is also strongly hy- 
perbolic). However, this point is of no consequence for 
the discussion that follows, so we will ignore it from now 
on. 

The evolution equation for F 1 given above is quite gen- 
eral, but it has the serious disadvantage that it involves 
quantities that are not tensors, such as T , mn and T* itself. 
In the next Section we will address this issue. 



B. Curvilinear coordinates 

The NOR formulation just described can in principle 
be used with any type of coordinates. However, when 
dealing with curvilinear coordinates, that is coordinates 
that are non-trivial even in flat space, one can easily find 
that the conformal connection functions T l are singular 
at some points and generally do not behave as a vec- 
tor would do. For example, in spherical coordinates the 
quantity r r turns out to be singular in flat space, while T e 
is non-zero even if we assume spherical symmetry. This 
in itself is not necessarily a major problem, as the equa- 
tions are quite general and are consistent in any set of 
coordinates. However, dealing with singular quantities 
numerically can be troublesome, and also dealing with 



non-tensor quantities makes it difficult to compare evo- 
lutions done with the same slicing conditions but differ- 
ent spatial coordinate systems. It would then seem like a 
good idea to replace the non-covariant quantities T l with 
a true vector. 

In order to do this we will start from the tensor A 1 jj c 
defined in equation (|3.ip above, and furthermore we will 
also assume that the background metric is the flat metric 
in the same curvilinear coordinates we are considering. 
Notice, in particular, that in Cartesian coordinates we 
have T l - k = 0, so that in that case A l jk and T l jk are 
identical. 

Just as we did before, we will again define the quan- 
tities A 1 as in (|3. 101) . We now want to calculate the 
evolution equation for A 1 . From the definition above we 
immediately find 

d t A l = dtP - f l mn d tl mn , (4.7) 

where we have used the fact that the flat background 
does not evolve. Using now (|4.5[) and the ADM evolution 
equation for 7^ one can easily find that 

d t A i -£ p A i = i mn d m d n p i + 1 mn £ p t i mn 

- [2K m - f m K] V m a 

- aV m [{2~OK lm -{l-Ol m K] 

+ 2aK™A i mn , (4.8) 

where the term £pt l mn must be calculated as if t l mn 
where a true tensor. In the previous equation A 1 is 
clearly a vector, and so is dt A 1 , but the right hand side is 
not manifestly covariant since it involves partial deriva- 
tives of the shift and terms containing T l mn . However, 
this can be easily fixed since one can show that, quite 
generally, 

7 m "V m V„/3 J = j mn d m d n (3 i + 1 mn £ r mn 

+ P^R'mnl , (4.9) 

with R l mn i the curvature tensor of the background. Since 
in our case the background is flat by construction, we can 
use the last result to rewrite the evolution equation for 
A 1 in the following way 

d t A* - £ & A l = 7 m ™V m V„/3 J 

- [2K im - f m K] V m a 

- aV m [(2-0K lm -(l-0l m K] 
+ 2aK mn A i mn . (4.10) 

The last equation is now manifestly covariant. 

In summary, in order to use the NOR formula- 
tion in curvilinear coordinates we need to express the 
3-dimcnsional Ricci tensor that appears in the evolution 
equations for the extrinsic curvature as (confront this 
with eq. flUQ)): 

Rab — - ^ l^^m^nlab + 7m( Q Vfc)A m + A m A( afc ) m 

+ 2A" m ( a A b )„ m + A mn a A mnb , (4-11) 
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with A a bc and A a denned in (|3~T1) and (|3TT0|) . pro- 
mote the A a to independent quantities, and evolve them 
through (j4~T0|) . 



V. THE BSSN FORMULATION 

A. Standard formulation 

The BSSN formulation is a reformulation of the ADM 
evolution equations, based on the work of Shibata and 
Nakamura [f| and Baumgarte and Shapiro [f|, that has 
proven to be particularly robust in the numerical evolu- 
tion of a large variety of spacetimes. This formulation is 
based on a conformal decomposition of the metric of the 
form 



7y = e 4< Hj 



(5.1) 



where the conformal factor is chosen in such a way that 
the determinant of the conformal metric is unity 7=1, 
which implies: 



= — In 7 
12 ' 



(5.2) 



From the definition above and the ADM evolution equa- 
tion for the spatial metric (12.71) . one can easily find the 
following evolution equation for </>: 



d t <f> 



(aK - d m p m ) + P m d m ct> 



The last equation can in fact be rewritten as 

1 



d t 4> - £p(f> 



6 



aK 



(5.3) 



(5.4) 



(5.5) 



where the Lie derivative of <f> is given by 

£p4> = P m d m $ + I d m (3 m . 

6 

Notice that, strictly speaking, (f> is not a true scalar 
density since its definition involves a logarithm, but 
ijj := — 7 1 / 12 is a well defined scalar density of weight 
1/6, so that the Lie derivative of </> is just £p(j) = £p' t P/' t P> 
which reduces to the expression above. 

The BSSN formulation also separates the extrinsic cur- 
vature into its trace K and its trace-free part 



Aij = Kij - - - , , K 



(5.6) 



We further make a conformal rescaling of the traccless 
extrinsic curvature of the form 



Aij — e — e ^ ( Kij — jijK 



(5.7) 



Just as we did in the case of the NOR formulation, the 
BSSN formulation also introduces three auxiliary vari- 
ables known as the conformal connection functions and 
defined as 



^ k f) k 



(5.8) 



where T l jk are the Christoffel symbols of the conformal 
metric, and where the second equality comes from the 
fact that the determinant 7 is equal to 1. 

The evolution equation for cj> was already found above, 
while those for 7^, K and A^ can be obtained directly 
from the standard ADM equations. The system of evo- 
lution equations then takes the form 



dtjij - £p% = -2aAij , 

d t <t> - £f3<f> = --aK , 
b 

d t A ij -£pA i j = e-^{-ViVja + aRij 



(5.9) 
(5.10) 



, TF 



Ana [ llJ (S-p)-2S lJ }} } 
a (KAij - 2A lk A k ^j , (5.11) 

i), K - £pK = -V' J n +n ( ,l,,_r' +-K 



+ 4ira(p + S) , 



(5.12) 



with V 2 := V m V m the spatial Laplacian operator associ- 
ated with the full physical metric, and where TF denotes 
the trace-free part of the expression inside the brack- 
ets. Notice also that indices of conformal quantities are 
assumed to be raised and lowered with the conformal 
metric. Here it is important to mention that the Hamil- 
tonian constraint has already been used in the evolution 
equation for K in order to eliminate the Ricci scalar. 

In the evolution equation for Aij above one needs to 
calculate the Ricci tensor associated with the physical 
metric, which can be separated into two contributions in 
the following way: 



R; 



(5.13) 



where Rij is the Ricci tensor associated with the confor- 
mal metric 7^ , which we write in terms of the T l as 



Rij = -fl mn d m d n % +l k(l d 3) T k + T k T (l0)k 



(* j)mn 1 ^ mnj j 



(5.14) 



(this is just the standard expression (|2.9[) for the con- 
formal metric), and where Rfj denotes additional terms 
that depend on derivatives of <j>: 

R% = -2ViV^ - 27y V fe V fc 

+ 4V^V^-47^V fc 0V fe , (5.15) 

with Vj the covariant derivative associated with the con- 
formal metric. 

Notice also that the evolution equations for Aij and K 
involve covariant derivatives of the lapse function with 
respect to the physical metric 7^ (i.e. covariant deriva- 
tives with no tilde). One must also be careful with the 
fact that in the evolution equations above we need to cal- 
culate Lie derivatives with respect to the shift vector /3 l 
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of tensor densities. In particular, -yy and Ay are tensor 
densities of weight —2/3. 

We are still missing an evolution equation for the T\ 
This equation can be obtained directly from the defini- 
tion, equation (|5.8I) . One finds: 

d t f i -£ r i = fk djdk/3 i + ^~ ijdjdkP k 

- 2(n,V ; ,l-' • ,i" ( V.n) . (5.16) 

In the above equation the Lie derivative of P should be 
calculated as if T l where a vector density of weight 2/3: 

£ p r = /Pdjf* - T J d 3 p l + \ rdjF . (5.17) 

Again, just as we did in the case of NOR, we will mod- 
ify the evolution equation for T* given above by adding 
to it a multiple of the momentum constraints. In order 
to do this, let us first rewrite the Hamiltonian and mo- 
mentum constraints in terms of the conformally rescaled 
quantities. One finds: 

H := \ \R+ \k 2 - AijA^j - 87rp = , (5.18) 

M' : VjA ij - - 8nf 

~ 87Tj l = . (5.19) 

where Vj now denotes covariant derivative with respect 
to the conformal metric. Notice also that the fact that 
the covariant metric has unit determinant implies that 
the term VjA J can be written as: 

VjA ij = d 3 A 11 + F jk A sk . (5.20) 

The Hamiltonian constraint above was in fact already 
used in order to eliminate the Ricci scalar from the evo- 
lution equation for the trace of the extrinsic curvature K 
above (equation (|5.12[) ]). 

Adding now a multiple of the momentum constraint to 
the evolution equation for T l , equation (|5.16[) above, we 
find: 

ftf - £ fi f 1 = f k 8 3 d k (3 l + 1 T ii d j d h p k - 2A ij d a 

- a(2-C)9 j l ii +«C(fM ii 

+ SA^drf - p ij djK - 87rj^ . (5.21) 

with £ an arbitrary parameter, and where j l := e 4 ^j\ 
The standard BSSN formulation usually takes £ = 2, 
which seems to be an optimal choice. 

Just as in the case of the NOR formulation, The BSSN 
formulation just described can be shown to be strongly 



hyperbolic for £ > 1/2 @. Standard BSSN with £ = 2 
has turned out to be particularly robust in practice, 
and leads to stable and well behaved numerical simu- 
lations. In conjunction with the Bona-Masso slicing con- 
dition (14. 6|) , and the so-called "Gamma driver" shift con- 
dition [ljj (see Section fV Cl below). it has allowed for the 
accurate simulation of the inspiral collision of black holes 
with different masses and spins [7-9, 18]. Today, most 
3-dimensional production numerical relativity codes use 
the BSSN formulation in one way or another, the notable 
exception being codes that use the "generalized harmonic 
formulation" (see e.g. [Hj]). 

B. Curvilinear coordinates 

When adapting the standard BSSN formulation to 
curvilinear coordinates we are faced with two problems. 
The first one is essentially the same problem that we had 
with the NOR formulation, namely that the quantities 
T l are not vectors (or more specifically vector densities, 
but we will come back to that point below). The sec- 
ond problem is the fact that in curvilinear coordinates 
the determinant of the flat metric is generally different 
from unity, so that asking for 7 = 1 is not a good idea. 
Consider, for example, flat space in spherical coordinates 
(r, 0, ip) for which the spatial metric is: 

ds 2 = dr 2 + r 2 dQ 2 , (5.22) 

with d£l 2 = d0 2 + sin 2 Odtp 2 the standard solid angle ele- 
ment. We then find that 7 = r 4 sin 2 9. 

In curvilinear coordinates it is in fact much better to 
ask for the determinant of the conformal metric to reduce 
to its value in flat space (see e.g. [l|). In order to avoid 
confusion between the conformal metric of the standard 
BSSN formulation and the one we will use here, from 
now on we will denote conformal quantities with a hat 
instead of a tilde. Also, for the conformal factor we will 
use x instead of (j), so that we will in fact have 

7y = e~ 4x ~fij , (5.23) 

and we will ask for j(t — 0) = 7, with 7 the determinant 
of the flat metric background in the same curvilinear co- 
ordinates. This change introduces two new features into 
the BSSN formulation. In the first place, in general we 
will find that 7 will not be constant in space so that we 
can no longer ignore its spatial derivatives. But more im- 
portantly, it is now not immediately clear how 7 should 
evolve in time. 

Following Brown |l) , one can suggest at least two "nat- 
ural" choices for the evolution of 7: 

1. dtj = 0. This is called a "Lagrangian" condition 
since the determinant of the conformal metric is 
constant along time lines. 

2. dtj— £/37 = 0. This is instead an "Eulerian" con- 
dition, since the determinant of the conformal met- 
ric is now constant along the normal lines (i.e. it 
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remains constant in time as seen by the Eulerian 
observers), so that it can in fact evolve along the 
time lines. 

Standard BSSN then corresponds to the Lagrangian 
case in Cartesian coordinates. Using now the fact that: 

£pl = P m d m l + 2jd m f3 m = 2jV m (3 m , (5.24) 

we can write in general for the evolution of 7: 

da = s (2 7 V m /3 m ) , (5.25) 

with: 

Lagrangian , 



1 Eulerian 



(5.26) 



On the other hand, since now 7 ^ 1, we find for the 
conformal factor \: 

1 



X = 12 ln ^ 7 ^ 



(5.27) 

We can now use this to find the evolution equation for x- 



1 f d t j da 

dtX = 777 ( ~ 

12 V 7 7 



= ^-UaK+^l-s^L) , (5.28) 
12 V 7 7 



which implies: 



6 6 
6 6 



(5.29) 



where we have used equation (j5.24j) above, and where 
£pX '■— £pl ' ll ~ £plll- In the above equation we have 
also introduced the shorthand a = (1 — s), so that a = 1 
now corresponds to a Lagrangian evolution and a = to 
an Eulerian evolution. 

There is an important point regarding the tensorial 
character of x that should be mentioned here. Notice 
that because of the new definition of x, equation f|5.27j) . 
we now have: 



„ 1 (£pn £pi 

Z 3 ^ — To 



12 V 7 7 



12 



V m /3 m - V m /3 
/3 m 9 m ln( 7 / 7 ) , 



so that finally 



£ P x = rdmx ■ 



(5.30) 



(5.31) 



In other words, the Lie derivative of x is now that of 
a scalar function with no density weight. We will see 



below that this will be the case for all dynamical quan- 
tities. This is another important difference between the 
standard BSSN formulation and the generalization we 
are introducing here, and it can be traced back to the 
fact that the definition (|5 . 2T[) of the conformal factor x 
now involves the ratio of two volume elements, so that x 
is a true scalar. 

The next step is to find the evolution equation for 7^ . 
Starting from the definition (|5.23|) above, and using the 
evolution equation for x (15.29[) . together with the ADM 
evolution equation for 7^ given by (I2.7[) . we now find 



dtlij ~ £p% 



-2a A, 



where A, is now defined as: 



1 



A tj : e- 4 * ( - JijK 



(5.32) 



(5.33) 



Again, with the definition of x above, both 7^ and Aij 
are true tensors and not tensor densities, and their Lie 
derivatives should be calculated accordingly. 

Similarly, the evolution equations for Aij and K be- 
come: 



d t Aij - £$Aij 



e- ix {-V l V J a + aR ij 
4ira (S - p) - 2Sij]} TF 
a (KAij - 2A ik A k j 
2 



-aA l3 V m f3 m , 



(5.34) 



d t K - £ f3 K = -V z a + a \ AijA 13 + -K 
+ 4:Tra(p + S) , 



1 



(5.35) 



Again, the Lie derivatives on the left-hand side are now 
those of proper tensors with no density weight. Notice 
that there is no term with crV m /3 m in the evolution equa- 
tion for K since it is a scalar. 

Just as before, the Ricci tensor that appears in the 
evolution equation for Ay is now separated into two con- 
tributions in the following way 



Rij — I!,, + Rfj 



(5.36) 



where Rij is the Ricci tensor associated with the con- 
formal metric 7,., , and where Rf- denotes the terms that 
depend on derivatives of x : 



R x 



= -2ViVjX - 2% V fc V fc x 

+ 4V lX V jX -47yV fc xV fc x, (5-37) 



with Vj the covariant derivative associated with the con- 
formal metric 7y . 

Following what we did in the case of the NOR formu- 
lation, we now want to write a fully covariant expression 
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for the conformal Ricci tensor R\j . In order to do so we so that the evolution equation for A 1 takes the final form 
will again introduce the quantities 

A a 6c := f \ c - f a bc , (5.38) 
A* := f nn A\ nn = r - 7 m "f \ nn . (5.39) 



d t A* - £ A t = 7 m "V m V„/3 l 



With these definitions the conformal Ricci tensor can be 
written as 



Rab 



1 



7 m "V m V„ 7 a 6 + 7 m(a V fc) A™ + A m A (ab)r 



2A mn (a A b)mn + A™"A 



mnb ; 



(5.40) 



The next step is to promote the A 1 to independent 
variables and find an evolution equation for them. In 
order to do this we must first find the evolution equation 
for the f \ Notice that, since now we have 7^1, the P 
now take the form 



L I 1 ran — u m I 



-f m a m ln7 



(5.41) 



Using equations (|5 -25[) and (|5 . 32[) we now find, after some 
algebra: 



(5.42) 



+ | [f m d m (v„/3") + 2f l VnjS" 

where £pT % is calculated as the Lie derivative of a vector: 

£pf l = p m d m t i - T m d m P . (5.43) 

The evolution equation for A 1 can now be obtained 
from the last equation using the fact that 



(5.44) 



where, just as we did in the case of the NOR formula- 
tion, we have assumed that the flat background does not 
evolve. One finds 



d t A l 



£pA i + r m d m d n p + r m £^ 1 

2V m ( aA m ) + 2aA mn A l mn 



V 



(V„j8 n ) + 2A l V„/3"] , (5.45) 



where again the Lie derivative of A 1 is that of a true vec- 
tor with no density weight, and the term £pT l m n should 
be calculated as the Lie derivative of a tensor. 

Just as it happened in the case of the NOR formula- 
tion, the right hand side of the last equations contains 
terms that involve partial derivatives of the shift, and 
also terms containing L l m „, so the expression is not ex- 
plicitly covariant. We can again fix this by using the fact 
that on a flat background the following relation holds: 



7 V m V„/5 =7 o m d n p +7 £pL 



(5.46) 



-0 

- 2V m (n,\ 
a 

+ 3 



■ 2aA" m A l 



V 1 (v„/3") + 2A 1 V„/3 n 



(5.47) 



which is now manifestly covariant. 

The last step is to add a multiple of the momentum 
constraints to the evolution equation for A* above. Doing 
that we finally find: 



7 m "V m V„/3 l - 2A m m a 
~im ,.2aA mn A i m „ 



d t A l - £ p A l 

V 2 (v„/3") + 2A J V n /3 



(5.48) 



where again ^ is an arbitrary constant that must be such 
that £ > 1/2 for the final system to be strongly hyper- 
bolic. 

We can now ask which would be the preferred choice 
for a in the above evolution equations, that is, should 
we take a Lagrangian or an Eulerian approach? Look- 
ing at the evolution equation for Xi equation (I5.29[) , one 
might first think that the simplest choice would be to 
take a = (s = 1), that is an Eulerian approach, since 
in that case the evolution equation simplifies. However, 
from the discussion above about the scalar character of x 
we see that if we choose a = 0, the evolution equation for 
X does not reduce to the standard BSSN evolution equa- 
tion for 4> given by equation (|5.4[) in the case of Cartesian 
coordinates (for which 7 = 1). This statement might 
seem somewhat puzzling since for a = equations (|5 .4[) 
and (|5.29p look identical. However, one must remember 
that x is a true scalar, while </> is a scalar density, so that 
their Lie derivatives are different. The same is true for 
the evolution equations of 7^, Aij and A*. 

It is in fact not difficult to convince oneself that if 
we want to recover the standard BSSN evolution equa- 
tions in the case of Cartesian coordinates we must choose 
a = 1, i.e. the Lagrangian approach, and simply remem- 
ber that all dynamical quantities are now true tensors 
with no density weight. The terms corresponding to the 
Lie derivatives of tensor densities in standard BSSN now 
appear explicitly on the right-hand side of the evolution 
equations through the terms proportional to V m /3 m . 



C. The Gamma driver shift condition 

The equations presented in the previous Section are 
completely general, and can be used with any gauge con- 
dition, both for the lapse and the shift. One we men- 
tioned the issue of hyperbolicity we specified a particular 
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slicing condition (the Bona-Masso condition), and the use 
of a shift known a priori, but this was just in order to 
make the discussion concrete. 

One particularly important shift condition that has 
turned out to be extremely robust in practice, and in the 
last few years has allowed for the stable and accurate sim- 
ulation of inspiraling black holes is the so-called "Gamma 
driver" shift condition [17| . This shift condition is partic- 
ularly well-adapted to the BSSN formulation, and comes 
in two versions (with several variations). The first possi- 
bility is the "parabolic" Gamma driver which takes the 
form: 

dtF = cxdtt 1 , (5.49) 

with ci some positive constant. The reason for the name 
"parabolic" is that, since the time derivative of T l that 
appears in the right-hand side involves second derivatives 
of the shift, the above equation results in a generalized 
heat-like equation for the shift components. In numerical 
simulations, the above condition has the same problem as 
any other parabolic equation, namely that for numerical 
stability the time step must be proportional to the square 
of the spatial grid spacing, resulting in a prohibitive use 
of computational resources in the case of high resolution 
simulations. 

The second version of the Gamma driver is the so- 
called "hyperbolic" Gamma driver which can be written 
in two alternative forms. The first form is simply 

d t p l = c 2 F , (5.50) 

while the second is 

d 2 p l = c 2 d t r . (5.5i) 

In both cases we end up with a generalized wave equa- 
tion for the shift, which justifies the name "hyperbolic" . 
The second version is usually preferred since it allows one 
to add a damping term that has been found to be very 
important in numerical simulations: 

d 2 t p l = c 2 d t r - ndtP 1 . (5.52) 

For reasons that we will not go into here, typical values 
of the parameters are c 2 = 3/4 and rj = 2/Madm (notice 
that rj has dimensions of inverse distance, so it is usually 
scaled with the total ADM mass of the spacetime). 

The main problem with the above shift conditions from 
the point of view of our present discussion is that, while 
they work well in Cartesian coordinates, they are seri- 
ously flawed in curvilinear coordinates since on the left- 
hand side we have a proper vector, while on the right- 
hand side we have contracted Christoffcl symbols. But 
we see that this problem can now be easily solved by 
choosing conditions of the form 

d t /3 l = Cl d t A l parabolic , (5.53) 
d 2 /3 l = c 2 d t k l - rjdt/3 1 hyperbolic . (5.54) 



When working in curvilinear coordinates, one must then 
use these modified Gamma driver conditions (or rather 
"Delta driver" conditions) in order to keep everything 
consistent. 



VI. THE CASE OF SPHERICAL SYMMETRY 

Having found the general form of the NOR and BSSN 
equations in curvilinear coordinates, we will now consider 
the special case of spherical symmetry. Here we will con- 
centrate on the case of the BSSN formulation, for the 
NOR formulation the analysis is entirely analogous (and 
in fact somewhat simpler). 

A. BSSN in spherical symmetry 

1. Main equations 

We start by writing the general form of the spatial 
metric in spherical symmetry as 

dl 2 = e 4x (a( ri t)dr 2 + r 2 b{r, t)dn 2 ) , (6.1) 

with a(r, t) and 6(r, t) positive metric functions, dfl 2 the 
solid angle element dQ 2 = d9 2 + sin 2 Odip 2 , and where x 
is the BSSN conformal factor introduced in Section fVBI 
above. Notice that with this notation the components 
of the conformal metric are j rr — a, jgg — r 2 b, and 
% v = (r sinfl) 2 b. 

The determinants of the physical and conformal metric 
take the form: 

7 = ab 2 (r 4 e 12x sin 2 6) , (6.2) 
7 = ab 2 (r 4 sin 2 6) . (6.3) 

The determinant of the flat metric in spherical coordi- 
nates can be easily found by setting a = b = 1 in the 
expression for 7 above: 

7 = r 4 sin 2 9 . (6.4) 

The condition that j(t = 0) = 7 now implies that ini- 
tially we must ask for ab 2 = 1. 

Notice in particular that for a Lagrangian evolution 
[a = 1) the metric components a and b are in fact not 
independent of each other. This is because in that case 
the determinant of the conformal metric 7 remains con- 
stant in time, so that the relation ab 2 = 1 will always 
hold. In the Eulerian case (a — 0) the quantity ab 2 does 
evolve, but its evolution is entirely controlled by the shift, 
and is independent of both the lapse and the extrinsic 
curvature. 

Let us now consider the shift vector. Since we are 
in spherical symmetry, the shift (as well as any other 
vector) will only have a radial component: f3 z = (f3 r , 0, 0). 
Since the different evolution equations in BSSN involve 
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the conformal divergence of the shift it is convenient at 
this point to calculate it. One finds, after some algebra: 



V m (3 m = d r /3 r +l3 r 
= d r (3 r + /3 r 



d r a d r b 2 
IkT + ~b~ + r 
( drjab 2 ) 2 
I, 2ab 2 + r 



(6.5) 



Consider next the auxiliary vector A*. Since this is 
a true vector (this is the whole idea), it will again only 
have a radial component: A 1 = (A r , 0, 0). One can show 
that this in indeed the case from the definition (|5.39[) . 
For the radial component we find 



a 



d r a 
~2a 



d r b 
~b~ 



2 / a 
ri'-b 



(6.6) 



One must remember, however, that in what follows A r 
will be promoted to an independent variable and the 
equation above will be considered a constraint. 

Let us now find the specific form of the evolution equa- 
tion for the conformal factor \ an d the components of the 
conformal metric a and b. From equation (|5.29j) we find 
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(6.7) 



where the divergence of the shift is given by (|6.5p above. 
For the conformal metric components we find 



d t a = P r d r a + 2ad r j3 r 
d t b = p r d r b + 2b^ 



aa V m /3 m - 2aaA a , 



ab V m /3 m - 2abA b 



r 3 

where we have introduced the quantities 



(6.9) 



a: 



a,, 



A» 



(6.10) 



The reason for using the mixed components of the 
traceless extrinsic curvature instead of the fully covariant 
ones is that in spherical symmetry such a choice simpli- 
fies considerably the evolution equations. In particular, 
the fact that the tensor Ay must be traceless implies that 



A a + 2A b = . 
Notice also that in fact one has 

K = TArr = l Tr A rr = A\ 



(6.11) 



3.12) 



and similarly for the angular component, so when we use 
mixed components of second-rank tensors the "confor- 
mal" and "physical" versions are identical. 

Consider next the evolution equation for the trace of 
the extrinsic curvature K . We find: 



d t K = /3 r d r K ~V 2 a + < 
+ ina {p + S a + 2S b ) 



A, 



2A\ 



(6.13) 



with S a and S b the mixed components of the stress tensor 

S a := S r r , S b := S e e , (6.14) 

and where the physical Laplacian of the lapse is given by 
1 



-AX 



[d 2 r a 

d r a d r b 
~2a b~ 



(6.15) 



The evolution equation for the traceless part of the 
extrinsic curvature is somewhat more complicated. Re- 
member first that we only need an evolution equation for 
A a , since the traceless condition implies A b = —A a /2. 
Rewriting equation (I5.34[) for the case of spherical sym- 
metry we find 



t A a = /3 r d r A a - [ V r V r a - \w 2 a 



" \R r r - 



aKA a — 167ra (S a — S b ) 



where V 2 a was given above and 



V r V r a 



1 



d 2 a — d r a 



d r a 
~2a 



(6.16) 



(6.17) 



and where the mixed radial component of the Ricci tensor 
and its trace R are given by 



R', 



1 



d 2 a 
~2a 

2 



ad r A r 



2 \ b 2 



3 

" 4 

d r a 
rb 



d r a 
a 



Id-", 



4 d 2 X - 2d rX 



rd r b 



d r a d r b 



R 



Ax 

d r b 
~b~ 

1 - 



d 2 a d 2 b 
~2^ + ~b~ 

2 



2 
rb 



b 

ad r A r 
d r b 



(6.18) 



d r a 
a 



rX 



b 
d r a 
~2a 



s(d 2 x + (drx) 2 ) 

d r b 2 N " 
b r . 



(6.19) 



It is interesting to notice that in equation (|6.16[) there 
is no contribution from the divergence of the shift (and 
hence a plays no role). The reason for this is that even 
though such terms do appear in the evolution equation 
for A rri once we raise the index to find the evolution 
equation for A a — A r r they cancel out. In fact, the shift 
contribution reduces to a pure advection term j3 r d r A a . 
This is one of the reasons why working with the mixed 
components is useful. 
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Finally, we need an evolution equation for A r . Writ- 
ing (|5.48p for the case of spherical symmetry we find 



behavior for small r: 



d t A r = 



f3 r 8 r A r - A r d r P r + -d 2 r p r + j dri- 
ll b \ r 



6 \a 



(A a d r a + ad r A a ) 



2a ( A a A r ~^(A a - A b ) 



at 
a 



d r A a - - d r K + 6A a d rX 



(A a - A b 



d r b 
~b~ 



8%j r 



(6.20) 



with j r the (physical) covariant component of the mo- 
mentum density, and where as before £ is an arbitrary 
parameter such that £ > 1/2, with preferred value £ = 2. 

Finally, it is also convenient to write the specific form 
of the Hamiltonian and momentum constraints. One 
finds 



H = R-(Al + 2Af) +-K- 16np = , (6.21) 
M r = d r A a - % d r K + 6A a d rX 



+ (A a -A b )[- + 



drb 
~b~ 



8wj r = 



(6.22) 



2. Regularization 



As has already been discussed in [20, l21| , unless spe- 
cial care is taken, in spherical symmetry the coordinate 
singularity at the origin can be a source of serious prob- 
lems caused by the lack of regularity of the geometric 
variables there. The problem arises because of the pres- 
ence of terms in the evolution equations that go as 1 jr 
near the origin. At the analytic level, for a regular space- 
time one can show that such terms cancel exactly at the 
origin, thus ensuring well-behaved solutions. However, 
this exact cancellation usually fails to hold for numeri- 
cal solutions. One then finds that the numerical solution 
becomes ill-behaved near r = 0. 

There arc in fact two different types of regularity con- 
ditions that the geometric variables must satisfy at r = 0. 
The first type of conditions are simply those imposed by 
the requirement that the different variables should have a 
well defined parity at the origin, and imply the following 



a i~ 


- a + O{r 2 ) , 


(6.23) 


/T - 


- 0{r) , 


(6.24) 


a <- 


- a + 0(r 2 ) , 


(6.25) 


b r- 


- b + O(r 2 ) , 


(6.26) 


A a r 


^ K + 0(r 2 ), 


(6.27) 


A b r 


- Al + 0{r 2 ) , 


(6.28) 




- 0(r) , 


(6.29) 



with {a , a , b°, A®} perhaps functions of time, but 
not of r. 

Notice, however, that the above parity conditions are 
not enough to guarantee regularity of the system of equa- 
tions described in the previous Section. Although most 
terms involving divisions by powers of r are indeed man- 
ifestly regular given the different parity conditions (be- 
cause they involve various derivatives of the geometric 
quantities), there are in fact two types of terms that 
would seem to remain ill-behaved at the origin. In par- 
ticular, the expression for A r (equation (16.61) ) involves 
the term (1 — a/b)/r, while the expressions for the ra- 
dial component of the Ricci tensor and its trace R 
(equations (|6 . 1 8[) and (16. 19[) ^ have terms of the form 
(1 — a/b)/r 2 , which apparently blow up at the origin. 
Similarly, in the momentum constraint (|6.22l) we have a 
term of the type (A a — A b )/r, which again would seem 
to be ill-behaved at the origin. 

The reason why these apparently ill-behaved terms 
turn out to be regular after all is a consequence of a sec- 
ond type of regularity conditions. These new conditions 
come from the fact that spacetime should be locally flat 
at the origin, and imply that for small r we must have 



- b ~ 0(r 2 ) , A a ~ A b ~ 0(r 2 ) 



so that 



A, 



A 



(6.30) 



3.31) 



It turns out that it is not trivial to implement numeri- 
cally both the parity regularity conditions and the local 
flatness regularity conditions at the same time. The rea- 
son for this is that at r — we now have three boundary 
conditions for just two variables: both the derivatives of 
a and b must vanish, plus a and b must be equal to each 
other (and the same thing must happen for A a and A b ). 

The regularization issue has already been discussed in 
some detail in several references (20l - [22l |. Here we will 
introduce a regularization procedure based on the one 
presented in [20|, but with one important modification. 
We will then start by introducing an auxiliary variable 
defined as 



A:=4 fl-£ 



(6.32) 



Notice first that in [20| the variable A was in fact defined 
defined as A := (l — a/b)/r. This is just a small difference 
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of no real consequence. Now, the local-flatness regularity 
conditions above imply that close to the origin we must 
have 



A~ A° + 0(r 2 ) 



(6.33) 



The main difference with the regularization procedure 
described in [2(| is that we will now also introduce a 
second auxiliary variable defined as 



A x := 



1 



(A a - A b ) 



3.34) 



Again, the local-flatness regularity conditions imply that 
close to the origin we will have 



A x ~ A° x + 0(r 2 ) 



3.35) 



Having introduced A and A\ we can rewrite all appar- 
ently ill-behaved terms in the BSSN equations in terms of 
these quantities so that the equations now look regular. 
In particular, the expression for A r becomes 



d r a 
~2a 



d r b 
~b~ 



2rX 



3.36) 



while i?£ and R take the form 
1 



R', 



Ax 



d 2 a 
~~2a 

2 



ad r A r 



3 
' 4 

1 / d r b\ 1 , d r a 

2 — ~2 Adra+ ^b 



d r a 
a 



■2X11 



rd r b 



+ Ad 2 X -2d rX 



d r a d r b 



R 



Ax 



d 2 a 
~2a 

2 



a 
d 2 b 



(6.37) 



- a a.A r 



d r a 
a 



1 / d r b \ 'i i a 

+ Ht) + Vb{"-b 

+ 4X + 8(d 2 rX +(drX) 2 ) 
' d r a d r b 2 
~2a b~ r 



d r b 



rX 



(6.38) 



Similarly, the momentum constraint now becomes 



M r = d r A a ~-d r K + 6A a d rX 



Ax 2r 



d r b 



- SttjV = 



3.39) 



Of course, at this point we haven't really fixed the 
problem, all we have actually done is to define some new 
variables as short-hands where we have hidden the ill- 
behaved terms. The way to solve the problem is to pro- 
mote A and A x to independent variables (with initial data 
given through their definitions) , and find evolution equa- 
tions for them that are manifestly regular. 



The evolution equation for A can be found directly 
from its definition and the evolution equations for a and 
b (equations (|6.8|) and (|6.9|) ). and turns out to be 



d t X = /3 r d r X + - 
r 



o \ r 



2aa , .„ ,„. 
— A x . (6.40) 



Notice that this equation is manifestly regular as long 
as f3 r ~ 0(r) for small r, and Ax itself remains regular. 
Notice also that the equation does not involve a, so it has 
the same form for Eulerian or Lagrangian evolutions. 

In order to find the evolution equation for A x it is 
important to notice first that the traceless condition 
A a + 2Ab = implies that A x and A a are in fact related 
through 



.4, 



3A a 
2r 2 



(6.41) 



This actually implies that for ^4^ to remain regular we 
must ask for A a ~ 0(r 2 ) near the origin (or in other 
words: A° a = A° b = 0). 

We can now use the evolution equation for A a , equa- 
tion (|6.16[) , to obtain the evolution equation for A\ . After 
a long algebra, in which one needs to take care of keeping 
together several terms that might be ill-behaved individ- 
ually (particularly terms of the form d r (F/r), with F 
some function that behaves as F ~ 0(r) near the ori- 
gin), one finds 



d t A x = P r d r Ax + 2A X 



fir 



1 



rae 4 * 

a 
rae 4 * 

a 
"ae 4 x 

ctA 



2d, 



d r a\ d r a f d r a d r b 

r J 2r \ a b 

9rX \ _ ^rX ( ^rO_ , ^rb 

r J r \ a b 



^d*x+ a -d r[ - 

2a r \ r 



i 2b rb ~ 
r V 1+ a~Y A 

r \ b J a 

-olKAx — 8iraSx , 



d r a ( 3 d. 



4 a 



•SOrX 

Ad rX 



d r b 
~~b~ 



(6.42) 



where Sx '■= (S a — Sb)/r 2 . The above equation is now 
manifestly regular in all its terms. 

Notice that at this point one can in fact just choose 
to ignore the evolution equation for A ai evolve Ax us- 
ing (|6.42[) . and later recover A a though the relation 
A a = 2r 2 Ax/3. This is in fact what we do in the nu- 
merical code used for the numerical examples of the next 
Section, since in practice it seems to reduce considerably 
the size of the numerical error (though evolving both A a 
and Ax independently also results in regular and stable 
evolutions) . 
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VII. SOME NUMERICAL EXAMPLES 

We have constructed a numerical code using the reg- 
ularized BSSN formulation in spherical symmetry de- 
scribed in the last Section. The code uses a method 
of lines algorithm, with either second order 3-step itera- 
tive Crank- Nicholson (ICN) or fourth order Runge-Kutta 
(RK4) as the time integrator, and second or fourth order 
centered differences in space. This code turns out to be 
very robust, stable, and well behaved at the origin. 

Here we will present some examples of numerical sim- 
ulations using this code that involve both pure gauge 
dynamics in vacuum situations, and simulations with a 
non-zero matter field. We will also consider simulations 
of a Schwarzschild black hole, which is a special case since 
it is in fact not regular at the origin. 



A. Pure gauge dynamics 

As a first example we will consider a pure gauge pulse 
propagating through the numerical spacetime. We con- 
sider initial data corresponding to Minkowski spacetime: 



x = o, 

a = b = 1 , 

A a = A b = K = 

A r = , 



(7.1) 
(7.2) 
(7.3) 
(7.4) 



which also imply A = A\ — 0. Non-trivial gauge dy- 
namics are obtained by choosing an initial lapse with a 
Gaussian profile of the form 
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FIG. 1: Evolution of the trace of the extrinsic curvature K. 
The different panels correspond to times t — 0,5, 10, 15. 

Hamiltonian constraint 
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1 + r 2 



-(r-ro) 2 /* 2 +e -(r+ro) 2 /* 2 



(7.5) 



with ao some initial amplitude, r$ the center of the Gaus- 
sian and a its width. Notice that we have multiplied 
the whole expression with r 2 /(l + r 2 ) and have in fact 
added two symmetric Gaussians (centered at r = r$ and 
r = — rO). This is done in order to make sure that the 
initial lapse is both an even function of r, and vanishes 
at r = 0. Having set up this initial lapse we evolve the 
system using harmonic slicing (with zero shift): 



d t a 



-a 2 K . 



(7.6) 



For the simulation shown below we have chosen the 
initial data parameters as ao = 0.01, ro = 5, a — 1, with 
grid parameters given by Ar = 0.1 and At = Ar/2. Dur- 
ing the simulation the initial pulse first separates in two 
smaller pulses propagating in opposite directions (due 
to the time-symmetry of the initial data). The inward 
moving pulse later implodes through the origin at t ~ 5 
and starts moving outward. Figure [1] shows a snapshot 
of the evolution of the extrinsic curvature K at times 
t = 0, 5, 10, 15. Notice that at t = 5 the value of K at the 
origin is too large for the chosen scale (it reaches a value 
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FIG. 2: Evolution of the Hamiltonian constraint. The differ- 
ent panels correspond to times t = 0, 5, 10, 15. 



of ~ 0.1), however the evolution always remains well- 
behaved and the value of K at the origin later returns to 
zero as the pulse moves outward. 

Figure [2] shows the evolution of the Hamiltonian con- 
straint for this simulation. Notice again that it remains 
well-behaved as the pulse goes through the origin. A 
small remnant of constraint violation that does not prop- 
agate away can also be clearly seen centered at the initial 
position of the pulse r ~ 5. 

Finally, we will consider the issue of convergence of 
the code. Figure [3] shows the root-mean-square (RMS) 
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RMS norm of Hamiltonian constraint 



are then given by 
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FIG. 3: RMS norm of the Hamiltonian constraint as a func- 
tion of time, for simulations at three different resolutions: 
Ar = 0.1, Ar = 0.05 and Ar = 0.025. The two highest reso- 
lutions have been rescaled by factors of 4 and 16 respectively. 



norm of the Hamiltonian constraint as a function of time, 
for simulations at three different resolutions, Ar = 0.1, 
Ar = 0.05 and Ar = 0.025 (keeping always the same 
ratio At/ Ar = 0.5). For the two highest resolutions, the 
norms have been rescaled by the corresponding factors 
expected for second order convergence, 4 and 16 respec- 
tively. The fact that all three lines coincide indicates that 
the code is indeed converging to second order. 



Scalar field 



The second set of simulations are somewhat more in- 
teresting since they evolve a non-zero matter field and 
therefore contain true dynamics. The matter field will 
correspond to a real scalar field $ with stress-energy ten- 
sor given by: 



p := nVT^ = - II 



* 2 



ae 



-P r »n v T, 



1 



$2 



S h := TS = -[U 



J/2\ 



(7.10) 
(7.11) 
(7.12) 

(7.13) 



The scalar field evolves through the simple wave equa- 
tion □ <& = 0, which in this case can be written as the 
following system of first order equations: 



d t § = (3 r d r ^ + aU, 

d t V = /3 r d r V + <F<9 r /3 r + d r (all) 

d t n = p r d r u + [ dr y 

'2 d r a d r b 
r ~2a~ + T 



(7.14) 
(7.15) 



2d rX 



+ 



ae 



4X 



d r a + aKU 



(7.16) 



For the initial data, we again choose an initial Gaussian 
profile for the scalar field of the form 



$ = 



$or 2 
1 + r 2 



e -(r-r ) 2 /° 2 + e -(r+r B ) 2 /v 2 



(7.17) 



with $o the initial amplitude, r the center of the Gaus- 
sian and a its width. We also assume time symmetry, so 
that at t — we have Kij = and II = 0, which implies 
that the momentum constraint is identically satisfied. 

Finally, we choose a conformally flat metric with 
a = b = 1, and solve for the conformal factor tp = e x us- 
ing the Hamiltonian constraint, which in this case takes 
the form 



dftp + -d r tp + 2mP b p = . 



(7.18) 



Notice that when we substitute the value of p given above 
(with n = 0), this equation reduces to 



(7.7) 



Let us now consider the case of spherical symmetry, 
and define the auxiliary first order quantities: 



n := n^d^ = - (dt& - fi r d r $>) , (7.8) 

a 



* := <9 r $ 



(7.9) 



The corresponding energy density p, momentum density 
j r , and stress tensor Sij that appear in the 3+1 equations 



5 t V + - d r ip + (tt* 2 ) tP = 



(7.19) 



This is clearly a linear equation for tp, which simplifies 
considerably its numerical solution (we can simply use 
centered spatial differences and invert the resulting tridi- 
agonal matrix directly). 

As gauge conditions we have chosen zero shift and 
1+log slicing, which is given by 



d t a 



-2aK . 



(7.20) 



For the specific simulation discussed here, we have 
taken as initial data parameters: $o = 0.04, r = 5, 
(7 = 1. We have chosen this initial data to be rather 
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FIG. 4: Evolution of the scalar field 3>. The different panels 
correspond to times t = 0, 5, 10, 15. 



strong, but not quite strong enough to collapse to a black 
hole. This is on purpose since we are interested in the 
regularity at the origin, and the collapse of the lapse as- 
sociated with the formation of a black hole makes this 
issue less relevant as everything just freezes close to the 
origin. In Section IVII CI below we will consider the case 
of a single Schwarzschild black hole. 

In Figure 0] we show snapshots of the evolution of the 
scalar field <&, for a numerical simulation using a grid 
spacing Ar = 0.025 and time step At = Ar/2. One can 
clearly see how the initial pulse separates into ingoing and 
outgoing pieces. The ingoing part then implodes through 
the origin and starts moving out, though significantly 
deformed. 

Figure \E\ shows the evolution of the central value of 
the lapse function a as a function of time (actually the 
value at r — Ar/2 since the origin itself is staggered). 
Notice how at t ~ 7 the central value of the lapse drops 
below 0.2, indicating a very strong gravitational field. 
However, the lapse later bounces and returns towards 1, 
and no black hole forms. The simulation remains well 
behaved throughout, and the origin remains regular. 

Finally, in Figure [6] we show again a plot of the RMS 
norm of the Hamiltonian constraint for three different 
resolutions, Ar = 0.05, Ar = 0.025 and Ar = 0.0125 
(due to the large dynamical range, the plot is now log- 
arithmic). Again, the higher resolution runs have been 
rescaled by factors of 4 and 16 respectively. The fact 
that the lines lie on top of each other shows that the 
code is converging to second order, as expected. One can 
notice that from t ~ 7 to t ~ 15 the convergence is less 
than perfect and the lines do not align precisely. This 
behaviour is a reflection of the fact that the gravitational 
field is very strong so that very high resolution is needed 
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FIG. 6: Logarithm of the RMS norm of the Hamiltonian con- 
straint as a function of time, for simulations at three different 
resolutions: Ar = 0.05, Ar = 0.025 and Ar = 0.0125. The 
two highest resolutions have been rescaled by factors of 4 and 
16 respectively. 



to adequately capture the situation. 



C. Schwarzschild black hole 

As our final example we will choose a Schwarzschild 
black hole, so again we are back in vacuum. The 
Schwarzschild solution is static in standard coordinates, 
but these coordinates are ill-behaved at the horizon. We 
will therefore use isotropic coordinates in which the ini- 
tial spatial metric takes the form 



dl 2 = ijj 4 (dr 2 + r 2 dn 2 ) , 
where the conformal factor ip is given by 
ip = 1 +M/2r , 



(7.21) 



(7.22) 
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and with M the mass of the black hole. As is well 
known, the Schwarzschild solution in isotropic coordi- 
nates has the topology of a wormhole (Einstein-Rosen 
bridge), with the throat located at r — M/2 (coincident 
with the horizon at t = 0), and with a coordinate sin- 
gularity at r = which corresponds to the compactifica- 
tion of the asymptotic flat region on the other side of the 
wormhole. 

The presence of the coordinate singularity at r — im- 
plies that our regularization procedure is now wrong since 
it was based on the idea of space being locally flat at the 
origin, which is clearly not the case here. The origin is 
not even part of space, which is why this type of initial 
data is known as puncture initial data (space in fact cor- 
responds to 5ft 3 minus the point at the origin, so it is a 
"punctured" 5ft 3 ). The main consequence of this is that 
using the regularization procedure now fails and the sim- 
ulations quickly crash. However, we have found that if we 
simply turn off the regularization, and run without intro- 
ducing the variables A and A\, we can have very stable 
and accurate simulations with the exception of the first 
few grid points closer to the origin where the code fails 
to converge (but we do find convergence away from these 
points as the plots below show). A more careful look 
at the data shows that, for the simulations described be- 
low, close to the origin we have A a ~ r, so that A\ ~ 1/r 
(confront equation (|6.4ip ). which explains why the reg- 
ularization procedure fails. A more detailed analysis of 
the behaviour at the origin for black hole simulations is 
clearly needed, but this is outside the scope of this paper. 

For all the simulations shown here we have chosen 
maximal slicing. This corresponds to the condition 
K = dtK = 0, which leads to the following equation for 
the lapse function: 

2 (2 d r a d r b \ 

-aae ix [K l3 K lj + 4n (p + S a + 2S b )] = 0, (7.23) 

where K t] K^ = A\ + 2A\ + K 2 /3. Notice that this is 
again a linear equation for a, which can be solved by 
direct matrix inversion. The reason for choosing max- 
imal slicing is to make sure that the lapse collapses at 
r = 0, which would not happen with 1+log slicing. This 
is because the origin is in fact an infinite proper distance 
away, so that any slicing condition with a finite speed of 
propagation would never change the value of the lapse 
there. || 

For the shift we have chosen a Gamma driver condition 
of the form discussed in Section IV CI In the particular 
case of spherical symmetry this condition reduces to: 

8 2 f3 r = ^ d t A r - 7} d t (3* . (7.24) 

Here we have already chosen the coefficient of the term 
d t A r equal to 3/4 in order to have an asymptotic gauge 
speed equal to 1. The condition above is solved in first 
order form by introducing the time derivative of the shift 



as an auxiliary quantity, so that we in fact solve the sys- 
tem: 

d t p r = B r , (7.25) 
d t B r = ~d t A r -r)B r . (7.26) 

Notice that this is the same shift condition (with minor 
variations) that is currently being used in most 3D codes 
that evolve black hole spacetimes. In the simulations 
shown below, the damping coefficient is always taken to 
be r) = 2. 

We still need to mention one final ingredient that 
goes into these simulations. Following [23, we have 
found that the simulations are better behaved if instead 
of evolving the singular conformal factor x directly, we 
evolve the quantity X := e~ 2x . [26] 

We are now ready to describe the numerical simula- 
tions. In all our simulations we have chosen the mass 
of the black hole to be M = 1 . We have chosen a grid 
spacing of Ar = 0.01 and time step of At = 0.005. We 
have also used 10,000 grid points in order to place the 
boundaries sufficiently far away so as not to have large 
errors from the boundaries affect the evolution. [27| 

Figures ITIfTOl show snapshots of the evolution of the 
lapse function a, the radial component of the shift vector 
ff , the radial metric component a and the conformal 
factor x at times £ = 0,5,10,15. Do notice that in order 
to better appreciate the plots, for the conformal factor 
X we have used a log plot, while for the shift we plot a 
smaller radial domain. 

The first thing to notice from these plots is the fact 
that the simulation is well behaved. The lapse collapses 
to at the origin (but with a non-zero derivative there as 
expected from the gauge conditions used, see e.g. 23}). 
The shift grows with time to counteract the slice stretch- 
ing effect, but becomes almost stationary at late times. 
While the conformal factor, though still singular at the 
origin also remains well behaved. 

In order to show that the simulations remain well be- 
haved for long times, and in fact reach an almost station- 
ary state, in Figure[TT]below we show the evolution of the 
maximum value of the radial metric component a and the 
radial shift vector /3 r up to t = 100. In both plots we can 
see that initially the maximum values grow rapidly, but 
this behaviour is later replaced by a very slow upward 
drift (this drift is well known from 3D simulations and is 
a consequence of the gauge conditions, particularly the 
damping term in the Gamma driver shift condition). 

Figure fTJ] shows the time evolution of the coordi- 
nate position of the apparent horizon and the appar- 
ent horizon mass (defined in terms of its area A a h as 
Mah = {Aah/l&ir) 1 / 2 ). One can notice how the radial 
position of the horizon drifts outward from r = 0.5 ini- 
tially to r ~ 1.1 at the end of the simulation. On the 
other hand, the horizon mass remains within 0.005% of 
unity throughout the entire simulation. 

This shows that the spherically symmetric BSSN code 
with maximal slicing and a Gamma driver shift condi- 
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FIG. 7: Evolution of the lapse function a. The different pan- 
els correspond to times t = 0, 5, 10, 15. 



FIG. 9: Evolution of the radial metric function a. The differ- 
ent panels correspond to times t = 0, 5, 10, 15. 
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FIG. 8: Evolution of the radial component of the shift vector 
/3 r . The different panels correspond to times t = 0,5, 10, 15. 



FIG. 10: Evolution of the log of the conformal factor \- The 
different panels correspond to times t = 0,5, 10, 15. 



tion can successfully and accurately simulate a black hole 
spacetime for a very long time without the need to excise 
the black hole interior. 



VIII. CONCLUSIONS 

Following Brown [lj, in this paper we have described 
how to modify standard hyperbolic formulations of the 
3+1 evolution equations of General Relativity in such 
a way that all auxiliary quantities are true tensors, 



thus allowing for these formulations to be used with 
curvilinear sets of coordinates. We have considered in 
particular both the Nagy-Ortiz-Reula (NOR) and the 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formula- 
tions, but the main ideas presented here can in principle 
be applied in general. 

The key idea has been that instead of using contracted 
Christoffcl symbols as auxiliary variables, one should use 
the difference between the physical Christoffel symbols 
and the Christoffel symbols associated with a background 
flat metric in the curvilinear coordinates, since the dif- 
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FIG. 11: Evolution of the maximum value of the radial metric 
component a (upper panel), and the radial component of the 
shift vector /3 r (lower panel). 
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FIG. 12: Coordinate position of the apparent horizon (up 
per panel), and apparent horizon mass M a h = (Aah/lftir) 1 ^' 
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ference of Christoffel symbols always behaves as a true 
tensor. 

Also, in the particular case of the BSSN formulation, 
one should not force the conformal volume element (the 
metric determinant) to be equal to 1, but rather to be 
equal to its initial value in the curvilinear coordinates. 
One important consequence of this choice is that all dy- 
namical geometric quantities now remain true tensors in- 
stead of tensor densities as in standard BSSN. One is 
also free to choose how this conformal volume element 
will evolve in time. Two "natural" choices present them- 
selves: the Lagrangian approach where the conformal vol- 
ume elements remain constant along time lines, and the 
Eulerian approach where they remain constant along the 
normal direction to the hypersurfaces. Standard BSSN 
can then be shown to be equivalent to the Lagrangian 
approach. 

Having developed the general formalism, we consid- 
ered as an example the particular case of BSSN in spher- 
ical symmetry, and studied in some detail the important 
problem of the regularity of the equations at the origin. 
For this we introduced extra auxiliary quantities that al- 
lowed us to impose the "local flatness" regularity condi- 
tion in a consistent way. Our regularization algorithm 
is similar to the one presented in [2(3], but it has been 
modified in a way that makes it more general and eas- 
ier to implement. It is important to mention that this 
regularization assumes that spacetime is regular at the 
origin and as such does not work for the case of black 
hole spacetimes where the origin is in fact a compactifi- 
cation of an asymtptic infinity on the other side of the 
Einstein-Rosen bridge. 

Finally, we presented a series of numerical simulations 
of our BSSN code in spherical symmetry, and we showed 
that the code was capable of evolving both regular space- 
times (with and without matter), as well as black hole 
spacetimes in a stable, accurate and robust way. 
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